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Abstract 

It has been recently argued through numerical work that rotating stars with a 
high degree of differential rotation are dynamically unstable against bar-mode de- 
formation, even for values of the ratio of rotational kinetic energy to gravitational 
potential energy as low as O(0. 01). This may have implications for gravitational 
wave astronomy in high-frequency sources such as core collapse supernovae. In this 
paper we present high-resolution simulations, performed with an adaptive mesh re- 
finement hydrodynamics code, of such low T/|W| bar-mode instability. The complex 
morphological features involved in the nonlinear dynamics of the instability are re- 
vealed in our simulations, which show that the excitation of Kelvin-Helmholtz-like 
fluid modes outside the corotation radius of the star leads to the saturation of the 
bar-mode deformation. While the overall trends reported in an earlier investigation 
are confirmed by our work, we also find that numerical resolution plays an impor- 
tant role during the long-term, nonlinear behaviour of the instability, which has 
implications on the dynamics of rotating stars and on the attainable amplitudes of 
the associated gravitational wave signals. 

Key words: gravitational waves, hydrodynamics, instabilities, stars: neutron stars: 
rotation 
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1 Introduction 



Neutron stars following a core collapse supernova are rotating at birth and 
can be subject to various nonaxisymmetric instabilities (see e.g. [1] for a re- 
view). Among those, if the rotation rate is high enough so that the ratio of 
rotational kinetic energy T to gravitational potential energy W, (3 = T/|W|, 
exceeds the critical value f}& ~ 0.27, inferred from studies with incompressible 
Maclaurin spheroids, the star is subject to a dynamical bar-mode (I — m — 2 



Preprint submitted to Elsevier 



1 February 2008 



/-mode) instability driven by hydrodynamics and gravity. Its study is highly 
motivated nowadays as such an instability bears important implications in 
the prospects of detection of gravitational radiation from newly-born rapidly 
rotating neutron stars. 

Simulations of the dynamical bar-mode instability are available in the litera- 
ture, both using simplified models based on equilibrium stellar configurations 
perturbed with suitable eigenfunctions [2,3,4,5], and more involved models for 
the core collapse scenario [6,7,8,9], and in either case both in Newtonian grav- 
ity and general relativity. Due to its superior simplicity the former approach 
has received much more attention, notwithstanding that the conclusions drawn 
from perturbed stellar models may not be straightforwardly extended to the 
collapse scenario. 

Newtonian simulations of triaxial instabilities following core collapse were first 
performed by [6]. These showed that the bar-mode instability sets in when 
f3 ^> 0.27 and when the progenitor rotates rapidly and highly differentially 
Such conditions are met when the (artificial) depletion of internal energy to 
trigger the collapse is large enough to produce a very compact core for which a 
significant spun-up can be achieved. More recently, three-dimensional simula- 
tions of the core collapse of rotating polytropes in general relativity have been 
performed by [7]. These authors studied the evolution of the bar- mode insta- 
bility starting with axisymmetric core collapse initial models which reached 
values of (5 ~ 0.27 during the infall phase. These simulations showed that the 
maximum value of (5 achieved during collapse and bounce depends strongly 
on the velocity profile, the total mass of the initial core, and on the equa- 
tion of state. In agreement with the findings from the Newtonian simula- 
tions of [6], the bar-mode instability sets in if the progenitor rotates rapidly 
(0.01 < (3 < 0.02) and has a high degree of differential rotation. In addition, 
the artificial depletion of pressure and internal energy to trigger the collapse, 
leading to a compact core which subsequently spins up, also plays a key role 
in general relativity for a noticeable growth of the bar-mode instability. 

Whether the requirements inferred from numerical simulations are at all met 
by the collapse progenitors remains unclear. As shown by [10] magnetic torques 
can spin down the core of the progenitor, which leads to slowly rotating neu- 
tron stars at birth (~ 10 — 15ms). The most recent, state-of-the-art compu- 
tations of the evolution of massive stars, which include angular momentum 
redistribution by magnetic torques and spin estimates of neutron stars at 
birth [11,12], lead to core collapse progenitors which do not seem to rotate 
fast enough to guarantee the unambiguous growth of the canonical bar-mode 
instability. Rapidly-rotating cores might be produced by an appropriate mix- 
ture of high progenitor mass (M > 25M Q ) and low metallicity (N. Stergioulas, 
private communication). In such case the progenitor could by-pass the Red 
Supergiant phase in which the differential rotation of the core produces a 
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magnetic field by dynamo action which couples the core to the outer layers 
of the star, transporting angular momentum outwards and spinning down the 
core. According to [13] about 1% of all stars with M > 1OM will produce 
rapidly-rotating cores. 

On the other hand, Newtonian simulations of the bar-mode instability from 
perturbed equilibrium models of rotating stars have shown that (3d ~ 0.27 
independent of the stiffness of the equation of state provided the star is not 
strongly differentially rotating. The relativistic simulations of [5] yielded a 
value of /3 ~ 0.24 — 0.25 for the onset of the instability, while the dynamics 
of the process closely resembles that found in Newtonian theory, i.e. unstable 
models with large enough /3 develop spiral arms following the formation of 
bars, ejecting mass and redistributing the angular momentum. As the degree 
of differential rotation becomes higher Newtonian simulations have also shown 
that Pd can be as low as 0.14 [14]. More recently [15,16] have reported that 
rotating stars with an extreme degree of differential rotation are dynamically 
unstable against bar-mode deformation even for values of (3 of 0(0.01). 

Given its recent discovery and its potential astrophysical implications for post- 
bounce core collapse dynamics and gravitational wave astronomy, we present 
in this paper high resolution simulations of such low T/|W| bar-mode in- 
stabilities. This work is further motivated in the light of the few numerical 
simulations available in the literature. Our main goal is to revisit the simula- 
tions by [15] on the low bar-mode instability, and particularly to check 
how sensitive the onset and development of the instability is to numerical 
issues such as grid resolution. To this aim we perform Newtonian hydrody- 
namical simulations of a subset of models analyzed by [15] using an adaptive 
mesh refinement (AMR) code [17] which allows us to perform such three di- 
mensional simulations with the highest resolution ever used. Our simulations 
reveal the complex morphological features involved in the nonlinear dynamics 
of the instability, where the excitation of Kelvin-Helmholtz-like fluid modes 
influences the saturation of the bar-mode deformation. We advance that while 
the overall trends found by [15] are confirmed by our work, the resolution 
employed in the simulations does play a key role for the long-term behaviour 
of the instability and for the nonlinear dynamics of rotating stars, which has 
implications on the attainable amplitudes of the associated gravitational wave 
signals. We note that we plan to upgrade the existing AMR code to account 
for the effects of magnetic fields in order to attempt the current study in a 
more realistic setup. The present work is a step towards that goal. 

The paper is organized as follows: Section 2 gives a brief overview of the 
equations to solve. Their solution is outlined in Section 3 which also contains 
the bare details of the AMR code. The results of the simulations are discussed 
in Section 4. Finally Section 5 presents our conclusions. 
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2 Mathematical framework 



The evolution of a self-gravitating ideal fluid in the Newtonian limit is de- 
scribed by the hydrodynamics equations and Poisson's equation: 



dp 
dt 



+ V • (pv) = 



dv . 1 _ , 

— + v • V )v = — Vp-V0 
at p 

dE 

— + V-[{E + p)v] = -pvV0 

V 2 = AnGp 



(1) 

(2) 

(3) 
(4) 

where x, v = ^ = (v x ,v y ,v z ), and <f>(t, x) are, respectively, the Eulerian 
coordinates, the velocity, and the Newtonian gravitational potential. The total 
energy density, E = pe + \pv 2 , is defined as the sum of the thermal energy, 
pe, where p is the mass density and e is the specific internal energy, and the 
kinetic energy (where v 2 = v 2 + v 2 +v 2 ). Pressure gradients and gravitational 
forces are the responsible for the evolution. An equation of state p = p(p, e) 
closes the system. We use an ideal gas equation of state p — (T — l)pe with 
T = 2. 

The hydrodynamics equations, Eqs. (1-3), can be rewritten in flux-conservative 
form: 

du gf(u) dg(u) dh(u) 

where u is the vector of unknowns (conserved variables): 

u = [p,pv x ,pv y ,pv z ,E] . (6) 

The three flux functions F a = {f , g, h} in the spatial directions x, y, z, respec- 
tively, are defined by 



f (u) = pv x , pv 2 x + p, pv x v y , pv x v z , (E + p)v a 



(7) 



pVy, pV x Vy, PV 2 + P, pVyV Z , (E + p)Vy 



(8) 



h(u) = \pv Z) pv x v z , pv y v z , pv 2 z +p,(E + p)v z 
and the source terms s are given by 



(9) 
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Table 1 

Overview of the initial models and results of the simulations. The rows report 
the name of the model, the ratio of equatorial-to-polar radii (r e /r p ), the degree of 
differential rotation (A), the ratio of kinetic to potential energy (T/|W|), the size of 
the computational grid (L) and the location of the corotation radius (r c ) for the two 
resolutions used: high (AMR H) and low (AMR L). In models R1H and R2H the 
corotation radius lies outside the star. The real (frequency) and imaginary (growth 
rate) parts of the bar-mode 02 are shown, for the low and high resolution simulation 
in comparison with the numerical results and linear analysis by [15]. Note that for 
model D3 no linear analysis results are available. 



Model 




Dl 


D2 


D3 


Rl 


R2 


r e /r p 




0.805 


0.605 


0.305 


0.305 


0.255 


A 




0.3 


0.3 


0.3 


1.0 


1.0 


T/\W\ 




0.039 


0.085 


0.149 


0.253 


0.275 


L/r e 




4.06 


3.73 


3.21 


4.25 


4.03 


r c /r e 


AMR L 


0.38 


0.47 


0.58 








AMR H 


0.36 


0.48 


0.56 






Re{a 2 )/n 


AMR L 


0.76 


0.58 


0.41 








AMR H 


0.81 


0.55 


0.43 




0.82 




Shibata 


0.80 


0.60 


0.45 


0.92 


0.75 




linear 


0.80 


0.58 




0.92 


0.75 


Im(cr 2 )/^o 


AMR L 


0.0042 


0.0154 


0.0200 








AMR H 


0.0089 


0.0190 


0.0240 


0.0005 


0.1960 




Shibata 


0.009-0.013 


0.019-0.021 


0.013 


<0.002 


0.23 




linear 


0.015 


0.021 




<0.002 


0.20 



d(p d(p d(f> d(f> d(p d(f) 

> P~a ' - P o-» 'PIT' P Vx x P v y7T ~ P Vz '^~ 

ox ay az ox ay az 



System (5) is a three-dimensional hyperbolic system of conservation laws with 
sources s(u). 
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3 Numerical approach 



For our study of the low T/\ W\ bar-mode instability we perform high- resolution 
simulations of rotating neutron stars using a Newtonian AMR hydrodynamics 
code called MASCLET [17]. The implementation of the AMR technique in the 
code follows the procedure developed by [18]. The hydrodynamics equations 
are solved using a high-resolution shock-capturing scheme based upon Roe's 
Riemann solver and second-order cell reconstruction procedures, while Pois- 
son's equation for the gravitational field is solved using multigrid techniques. 
The accuracy and performance of the MASCLET code has been assessed in a 
number of tests [17]. We note that the code was originally designed for cos- 
mological applications, and here it is applied to simulations of self-gravitating 
stellar objects for the first time. 

The simulations are performed with two different grid resolutions. The low 
resolution grid consists of a box of size L with 128 3 zones, yielding a fixed 
resolution of L/128. We note that the effective resolution of our coarse grid 
is comparable to that used by [15]. Correspondingly, the high resolution grid 
consists of a base coarse grid of 128 3 cells, and one level of refinement composed 
of patches with maximum size of 64 3 cells (32 3 coarse cells). This yields a grid 
resolution on the finest grid of L/256. This resolution is enough to resolve the 
structures simulated, and hence no deeper refinement levels are needed. The 
patches are dynamically allocated covering those regions of the star where the 
highest resolution is required (highest densities). Typically only one patch is 
needed for spheroidal models, and 4-8 in models with toroidal topology. The 
use of AMR techniques in our high resolution simulations, allows us to save 
about a factor 4 in CPU time and memory with respect to a unigrid simulation 
with 256 3 cells. No symmetries are imposed in the simulations. To the best of 
our knowledge, in the investigations of the bar-mode instability performed by 
previous groups, grid resolutions as high as the ones we use here were never 
employed. 

As customary in grid-based codes [19,20] the vacuum surrounding the star 
is filled with a tenuous numerical atmosphere with density p/p max ~ 1CT 12 
and zero velocities, p max being the maximun density. Every grid cell with 
p/pmax < 10~ 6 is reset to the atmosphere values. A correct treatment of the 
atmosphere is essential for an accurate description of the stellar dynamics and 
correct computation of the growth rates of unstable modes. We have checked 
that values for the atmosphere higher than those we chose or a free evolution of 
the atmosphere altogether, lead to remarkable changes in the mode behaviour, 
growth rates, and frequencies. We have also checked that lower values for the 
atmosphere do not produce those changes, which ensures that our evolutions 
are not affected by the atmosphere values used in the simulations. 
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4 Results 

4-1 Initial data 

Differentially rotating stellar models in equilibrium are built according to the 
method of [21], and used as initial data for the AMR evolution code. The stars 
obey a polytropic equation of state P = Kp v with index T = 2. As [15] the 
profile of the angular velocity fl is given by 

(ro/r e ) 2 + A 2 ' 

where r e is the equatorial radius of the star, Qq is the central angular ve- 
locity, w is the distance to the rotation axis, and A parametrizes the degree 
of differential rotation, from A <C 1 for highly differentially rotating stars to 
A — > 00 for rigidly rotating stars. For comparison purposes these parameters 
are chosen as in some of the models of [15], and are summarized in Table 1. 
Models labelled D rotate with a high degree of differential rotation, as A = 0.3, 
and may therefore be subject to the low T/|W| bar-mode instability. We also 
consider models almost rigidly rotating, labelled R, prone to experience the 
"classical" bar-mode instability. Labels L and H in the models refer to low 
and high resolution respectively. 

Following [15] we perturb the initial density profile according to 

p-p' ' ( l + s ^f)' < 12 > 

the perturbation of the pressure given by the equation of state accordingly. 
A perturbation amplitude 5 = 0.1 is used in all our simulations. As we show 
below this form of the perturbation excites the I — m — 2 bar-mode. In 
addition, grid discretization can leak small amounts of energy to all other 
possible modes, which could in principle grow provided they were unstable 
and the simulations were carried on for sufficiently long times. 



4-2 Stability analysis 

To compare with [15] we calculate the distortion parameters i] + and i] x (and 
V = {rf + + rf x ) 1 / 2 ) defined as 
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(11) 
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Fig. 1. Power spectra of A m from m = 1 to m = 8 for model D3H. 
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where = x,y,z) is the mass-quadrupole moment 



(13) 



hj = J dx 3 px l x j . 



(14) 



For the study of the growth rate and interaction of the different angular modes 
within the star is useful to calculate the global quantity 



J rfx 3 p(x) 



-vmip 



(15) 



and Am = A m /A . We follow the time evolution of modes with m ranging 
from 1 to 8. Since our initial equilibrium models are axisymmetric and have 
equatorial plane symmetry, all A m are zero initially but once perturbed all 
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Fig. 2. Evolution of i] for models R1H (upper panel) and R2H (lower panel). Expo- 
nential fits to the peaks in the growing phase are overplotted as solid lines. 

initial models exhibit a dominant m = 2 component. Assuming that the modes 
behave as e - l (°>"-*- m ¥>) ; real part of a m can be obtained by Fourier trans- 
forming A m - In particular Re(<7 2 ), the bar-mode frequency, can be extracted 
from either A<i or r\ as both represent the same mode. This is the dominant 
mode in all our simulations and its frequency and growth rate are given in 
Table 1. The latter corresponds to the imaginary part of <72, which is calcu- 
lated fitting an exponential to the peak values of r\ in the growing phase of 
the evolution until the modes saturate. Other modes are also identified in the 
simulations for values of A m with lower amplitudes. We have checked that 
these modes are harmonics of the I = m = 2 mode so that they follow to good 
accuracy the relation a m = ma p , a p being the pattern frequency, calculated 
as ij p = (T2/2. This is shown for model D3H in Fig. 1 which displays the spec- 
trum of Am from m = 1 to m = 8 (in arbitrary units). The vertical dashed 
lines in this figure indicate the location of the integer multiples of the pattern 
frequency cr p , their values indicated on the axis at the top of the figure. Each 
spectrum for each mode is normalized to its own maximum for plotting pur- 
poses. Note that the lower the mode amplitude the noisier the spectrum and 
the less accurate the relation a m = ma p . 

For the models of our sample subject to the "clasical" bar-mode deformation 
(R1H and R2H), our simulations yield a value of /3 between 0.253 and 0.275, 
in good agreement with the critical value for the onset of the dynamical bar- 
mode instability. Model R1H is stable and model R2H is unstable. The growth 
rates and frequencies reported in Table 1 agree with those of [15]. Note that 
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R2H 




Fig. 3. Evolution of r\ for models Dl (upper panel), D2 (central panel) and D3 (lower 
panel). Dashed lines correspond to low resolution and solid lines to high resolution. 
Exponential fits to the peaks in the growing phase are overplotted as solid lines. 

for model R1H, which is stable, the frequency for the m = 2 mode cannot be 
computed. The time evolution of rj for these two models is displayed in Fig. 2. 
For the unstable model R2H, our simulations show the formation of a bar 
which saturates for values of rj + and i] x close to 1, i. e. in the full nonlinear 
regime. 

Fig. 3 shows the time evolution of 77 for models D in our sample, prone to suffer 
the low bar- mode instability. Solid lines correspond to high resolution 

simulations and dashed lines to low resolution. For all three models the pattern 
frequencies a p are such that there exists a corotation radius inside the star, 
i.e. a radius at which the bar-mode rotates with the same angular velocity as 
the fluid. The location of the corotation radius for all models of our sample 
is reported in Table 1. As recently discussed by [22] the existence of such 
corotation radius is a potential requirement for the ocurrence of the instability. 
As becomes clear from Fig. 3, all models are unstable but grid resolution has 
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an important effect on the saturation of the instability once the nonlinear 
phase has been reached, as well as in the long-term dynamics of the stars. 

In the linear phase of models D1H and D2H, the growth rates and frequencies 
agree with the results of [15] in both, the numerical simulations and the linear 
analysis (see Table 1). In the linear phase of model D3H, our frequencies are 
similar to the numerical results of [15], although our growth rates are about 
a factor two larger. We emphasize that no results are reported in the linear 
analysis for this model in the work of [15], and therefore this discrepancy can 
be an effect of the resolution used or of the characteristics of each numerical 
code. Increasing resolution leads to similar results in the frequencies but to 
higher growth rates. 

In the nonlinear phase, models Dl and D3 behave similarly for the two resolu- 
tions used (see Fig. 3), and also similarly to the results by [15] (compare with 
Fig. 3 of that paper). For model D2 we observe a radical change of behavior in 
the nonlinear phase of the mode evolution depending on the grid resolution. 
This has implications on the long-term dynamics of the star and, in particu- 
lar, on the attainable amplitudes of the gravitational radiation emitted, as we 
discuss below. 

It is worth mentioning the possibility that the unstable mode at the start of 
model D2H might excite some other mode in the corotation band, which could 
not otherwise be excited for lower grid resolution. As discussed by [23,24] in 
their study of differentially rotating shells, there are many zero-step modes in 
the band, so that the whole continuous spectrum could potentially be excited. 
In such case these modes would have very slow power-law growth. 

For all our models we have checked mass conservation along the evolution. The 
worst results are obtained for model D3H, for which mass is conserved within 
2.5% error when the instability saturates. At the end of the simulation (after 
48 orbital periods and 25000 iterations in the coarsest grid) the error has grown 
to only 6%. For all other models mass conservation is even more accurate. Note 
that these errors are within the round-off error of the code, and it is not related 
to the conservation properties of the numerical scheme itself. For a regular grid 
with 128 3 cells and a simulation employing 25000 iterations, the accumulated 
round-off error (binomial distribution) using single-precision arithmetics, is 
about V128 3 x 25000 x 10~ 8 = 0.0023 = 0.23%. Correspondingly, for a 256 3 
grid (with twice the number of iterations for the simulation) the error is about 
0.9%. Taking into account that this error affects the nonlinear evolution of the 
system, it is not surprising to have an error at the level of a few percent by 
the end of our high resolution simulations, for all conserved quantities. 

Figure 4 shows the evolution of A m for model D3 and for m ranging from 
1 to 8 for our two resolutions. According to this figure, the only two modes 
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Fig. 4. Evolution of \A m \ for model D3 with low resolution (top) and high resolution 
(bottom). The m = 2 mode is represented with thick solid line, m = 4 with thin 
solid line, m = 6 with dashed line, m = 8 with dot-dashed line, and all other odd 
m with dotted lines. 

relevant for the dynamics of the star are m = 2 and m — 4. All other modes 
have smaller amplitudes and play no role in the dynamics. Note that for odd 
m modes, the value of the integrated quantity A m , if close to zero, is extremely 
sensitive to very small numerical asymmetries, which are induced by the patch 
creation scheme of our AMR code. This explains the resolution differences in 
the initial values for odd m modes in Fig. 4 (at t — they start off at 10~ 8 
level for the low resolution simulation), although they saturate at the same 
value irrespective of the resolution. 

An important diagnosis for the accuracy of the results is the location of the 
center of mass during an evolution. The round-off error of the numerical code 
imposes controlled errors in mass and linear momentum, which results in tiny 
displacements of the center of mass. However small (one numerical cell in 
our runs) this unphysical displacement may hinder the correct analysis of the 
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Fig. 5. Effects of the artificial displacement of the center of mass (of only one 
numerical cell) on the time evolution of |^4i| for model D2H. The thin solid line 
shows a fictitiuos evolution resulting from the numerical artifact originated by the 
center of mass displacement. 

mode growth rates. For this reason all integrated quantities shown in Fig. 4 
are computed after correcting for the displacement of the center of mass, 
x ncw = x old — x CM , in a post-processing stage of the data analysis. Were this 
not done, a one-armed m = 1 mode would grow much faster than it should to 
bring up fictitious features in the plots. This is shown for model D2H in Fig. 5. 
The thick solid line in this figure corresponds to the evolution of the m = 1 
mode taking into account the correction for the center of mass displacement, 
while the thin solid line is the corresponding evolution of this mode without 
the correction. 



4-3 Gravitational waves 



The growth and saturation of the instability is also imprinted on the gravita- 
tional waves emitted. The gravitational waveforms h + and h x for models Dl, 
D2, and D3, computed using the standard quadrupole formula, are shown in 
Fig. 6. For a source of mass M located at a distance R those waveforms can 
be calculated from the dimensionless waveform amplitudes a + and a x as 



h sin 2 6M 2 

h+, x =a +iX — — — , (16) 

using G = c = 1 units. The resulting chirp-like signal in all the models, partic- 
ularly apparent for model D2L, indicates the presence of a bipolar distribution 
of mass within the star (see Sec. 4.4). 
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Fig. 6. Gravitational waves for models Dl to D3 extracted using the standard 
quadrupole formula. Thick (thin) solid lines correspond to low (high) resolution. 
Only the dimensionless waveform amplitude a+ is plotted. 

As mentioned before, the effects of grid resolution on the evolution of the 
nonlinear phase of the bar-mode are imprinted on the gravitational waveforms. 
Thick solid lines in Fig. 6 are the waveforms which correspond to the low- 
resolution models, and thin solid lines to the high-resolution counterparts. 
The evolution of r\ for model D3, displayed in Fig. 3, shows little deviations 
with grid resolution, and this translates into very similar gravitational wave 
patterns (bottom panel of Fig. 6), the differences becoming more noticeable in 
the nonlinear phase following saturation (Qot > 75). For model Dl (top panel), 
the differences also become more apparent at later times during the evolution, 
in good agreement with the dissimilar behaviour of the matter dynamics in 
this model, as encoded in the evolution of rj in Fig. 3. As happens for model 
D3 the first few cycles of the gravitational waveform, when the mode is still 
in the linear phase, are accurately captured for both resolutions. 

The major dependence of the waveform on the grid resolution is found for 
model D2. Again, the linear phase for the growth of the bar deformation is 
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Fig. 7. Snapshots of the density, vorticity, and specific angular momentum, for 
model D3H, at three representative instants of the evolution. All snapshots show 
slices of the stars in the equatorial plane. Quantities are normalized as follows: 
p/pmLc, r e w v /vi , and l v /(r e Vs ), where vi°^ is the initial velocity at the surface 
of the star. 



accurately captured irrespective of the resolution (and agrees with the per- 
turbative results of [15]). This is signalled in the perfect overlapping of both 
gravitational waveforms during the first three cycles (see the middle panel of 
Fig. 6). However, the different nonlinear dynamics of the bar-mode deforma- 
tion for this model, shown in the middle panel of Fig. 3, is severely imprinted 
on the gravitational waveform. Model D2H emits gravitational waves which 
have roughly one order of magnitude smaller amplitude than those computed 
for the corresponding low resolution model. 
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4-4 Morphology 



We next describe the morphological features encountered during the evolution 
of some representative models. Fig. 7 shows three snaphsots of the evolution 
of model D3H for the density (top), the azimuthal component of the vorticity, 
= (V x v) 9 (middle), and the specific angular momentum, I — r x v 
(bottom). From left to right the snapshots correspond to the initial time (flat = 
0), a time when the bar-mode instability is growing (Qot = 33.6), and the time 
when the instability saturates (f2 t = 70.4). Only the equatorial plane of the 
stars is shown in all these plots. Animations of all simulations performed are 
available at www . uv . es/ ~cerdupa/bars/. We note that our AMR code is able 
to dynamically place patches (e. g. between 4 and 8 in the D3H model) and 
evolves the system with continuous matching between patches, as exemplified 
in Fig. 7. 

The evolution of model D3H shows that as the m = 2 mode grows the star 
develops an ellipsoidal shape which remains spinning beyond saturation. Since 
the low j3 m = 2 mode saturates at lower values (77 ~ 0.1) than the classical 
bar-mode instability (77 ~ 1), no clear bars are visible in the density plot. At 
late times (Q t > 100) a "boxy" structure becomes apparent as the m = 4 
mode has grown to almost similar amplitude as the m = 2 mode (see anima- 
tions and Fig. 4). No other global features can be seen, consistent with the 
fact that \A m \ <C 1 for all modes other than m = 2 and 4. The vorticity plot 
shows that the m = 2 mode at Q t = 33.6 adopts the form of a two-armed 
spiral winding up around the central parts of the star. As the mode begins to 
saturate (Q t = 70.4) the spirals break apart into the outer layers in a turbu- 
lent flow reminiscent of the (shear) Kelvin-Helmholtz instability, and shock as 
they reach the atmosphere. These trends are also visible in the specific angular 
momentum plot. 

The presence of a corotation radius, at r/r e = 0.56 for model D3H, seems 
to play a role in the growth and saturation of the instability, in agreement 
with the recent findings of [25]. As the bar-mode grows, pressure waves carry 
angular momentum outside the corotation radius, which is deposited in the 
outer layers of the star. This excites Kelvin-Helmholtz-like instabilities in the 
fluid that break the mode outside the corotation radius. When this happens the 
m = 2 instability stops growing and no more angular momentum is extracted. 
Figure 8 shows late-time snapshots of the equatorial plane distribution of the 
density perturbation, i.e. (p — p^)/ p£^ x , f° r models D2H and D3H. The times 
are chosen well inside the nonlinear and saturation phase of the instability. 
This figure helps to interpret the mode dynamics and its saturation along the 
lines mentioned before: During the evolution the density perturbations are 
shed in waves from the center towards the outer layers of the star. At late 
times, when the instability saturates, such shedding stops, and the density 
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Qot= 150.6 100.4 




Fig. 8. Snapshots of the density perturbation at the equatorial plane for models D2H 
and D3H. The white solid curves indicate the location of the corotation radius. The 
white dashed boxes indicate the location of the patches for model D3H. 



perturbation reaches the largest values outside the corotation radius (depicted 
with white solid lines in Fig. 8), for either model. 

We note in passing that the corotation radius in all our high resolution models 
lies well inside the outer boundary of the finest box set up by the AMR 
refinement pattern, (see, e.g. the white dashed boxes depicted in the right 
panel of Fig. 8 indicating the location of the AMR patches for model D3H) 
This rules out the possibility of a numerical artifact resulting from the patch 
creation scheme of our AMR code being the cause for the different long-term 
evolution between low and high resolution models, particularly noticeable for 
model D2 in Fig. 3. 

Finally, Fig. 9 shows a comparison between models D2L and D2H at Qot = 101 
(i.e. well within the nonlinear phase), to highlight the effects of the numer- 
ical resolution on the morphology. From top to bottom this panel shows a 
schlieren plot (|Vlogp|) , w^, and I. The resolution differences in the evolu- 
tion of model D2 become apparent from this figure. In particular, the "boxy" 
structure becomes much more clearly visible in the low resolution simulation 
(D2L), indicating an excessive growth rate of the m = 4 mode. The presence 
of pressure waves is emphasized in the schlieren plot, very accurately captured 
in model D2H. Those waves, once the flow is driven to turbulence past the 
corotation radius, redistribute the angular momentum in the outer layers of 
model D2L in a much more pronounced way than for model D2H. 
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Fig. 9. Resolution comparison between models D2L and D2H once the instability 
has saturated. Only slices of the stars in the equatorial plane are shown. 

5 Summary and outlook 



We have presented AMR high- resolution simulations of the low T/|W| bar- 
mode instability of extremely differentially rotating neutron stars. Our main 
motivation has been to revisit the simulations by [15] on such instability, 
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assessing how sensitive the onset and development of the instability is to 
numerical issues such as grid resolution. We have addressed the importance 
of a correct treatment of delicate numerical aspects which may spoil three- 
dimensional simulations in (Cartesian) grid-based codes, always hampered by 
insufficient resolution, namely the handling of the low-density atmosphere sur- 
rounding the star, the correction for the center of mass displacement, and the 
mass and momentum conservation properties of the numerical scheme. Our 
simulations have revealed the complex morphological features involved in the 
nonlinear dynamics of the instability. We have found that in the nonlinear 
phase of the evolution, the excitation of Kelvin-Helmholtz-like fluid modes 
outside the corotation radii of the stellar models leads to the saturation of the 
bar-mode deformation. While the overall trends reported in the investigation 
of [15] are confirmed by our work, the resolution used to perform the simu- 
lations may play a key role on the long-term behaviour of the instability and 
on the nonlinear dynamics of rotating stars, which has only become apparent 
for some specific models of our sample (namely model D2). This, in turn, has 
implications on the attainable amplitudes of the associated gravitational wave 
signals. 

The work reported in this paper is a first step in our ongoing efforts of study- 
ing the dynamical bar-mode instability within the magnetized core collapse 
scenario. 
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